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Abstract 

A new algebraic cubature formula of degree 2n + 1 for the product 
Chebyshev measure in the d-cube with w n d /2 d ~ 1 nodes is established. 
The new formula is then applied to polynomial hyperinterpolation of de- 
gree n in three variables, in which coefficients of the product Chebyshev 
orthonormal basis are computed by a fast algorithm based on the 3- 
dimensional FFT. Moreover, integration of the hyperinterpolant provides 
a new Clenshaw-Curtis type cubature formula in the 3-cube. 



1 Introduction. 

A cubature formula with high accuracy is an important tool for numerical com- 
putation and has various applications. One of the applications is to construct 
polynomial hyperinterpolation, introduced by Sloan |17J . which is an approx- 
imation process constructed by applying the cubature formula on the Fourier 
coefficients of the orthogonal projection operator. 

A cubature formula of degree In — 1 with N nodes with respect to a measure 
da supported on a set takes the form 

[ p(x)dfi= J2 w tP® forall pe^ft), (1) 

where {w^}, called weights, are (positive) numbers, X n is a set of points, called 
nodes, 

£:=(&,&,...,&) ex„cfi (2) 

with card(A„) = N, and 11^ denoted the subspace of d-variate polynomials of 
total degree < m restricted to f2. For a cubature formula of degree 2n — 1 to 
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exist, it is necessary that 



N := card(X„) > dim(n^(0)) = ^+ ^ = ^(1 + o(l)). (3) 

There are improved lower bounds of the same order in terms of n. A challenging 
problem is to construct cubature formulae with fewer nodes, that is, with the 
number of nodes N close to the lower bound. 

In this paper we consider the case that the measure is given by the product 
Chebyshev weight function 

dp = w d (x) dx, w d (x) ■= -1 rr -L= (4) 

on the cube f2 := [—1, l] d . For d = 1, the Gaussian quadrature formula of degree 
2n — 1 needs merely N = n points. Our main result is a new family of cubature 
formulae that uses N w n d /2 d ~ 1 many nodes. For d — 2 these formulae are 
known to have minimal number of nodes. For d > 3 they are still far from 
the lower bound, but they appear to be the best ones that are known at this 
moment. We refer to Section 2 for further discussions. We present numerical 
tests on these cubature formulae in three variables and also apply them to 
constructing polynomial hyperinterpolation operator in three variables. 

For every function / £ C(£!) the /i-orthogonal projection of / on 11^ (f2) is 



S n f(x)= ^2 a a p a (x), a a := / f(x)p a (x)dn , (5) 

\a\<n Jn 

where x — (xi, X2, ■ ■ ■ , Xd) is a d-dimensional point, a is a <i-index of length |a| 

a= («!,..., a d ) eN d , \a\ :=a 1 + ... + a d , (6) 

and the set of polynomials {p a , < |a| < n} is any /x-orthonormal basis of 
II„(f2) with p a of total degree \a\ (concerning the theory of multivariate orthog- 
onal polynomials, we refer the reader to the monograph |9J). Clearly, S n p = p 
for every p £ H d (fl). Given a cubature formula ([lj of degree < 2n, we obtain 
from ([5| the polynomial approximation of degree n by the discretized Fourier 
coefficients {c a } 

f(x) W C,J(x) := ^ C aPa{x) , C a 1= ^ W^fiOPctiO , (7) 
|o|<n (,eX n 

where c a — a a and thus C n p = S„p = p for every p 6 n„(f2). This is the 
hyperinterpolation operator. It satisfies the basic estimate: for every / £ C(f2), 



||/-£ n /||^ (n) <2VM(n)^(/)^0, n^oo, (8) 

where E n (f) := inf {||/ — pW^ , p £ Il^(fi)}, so that it converges in mean. The 
convergence rate can be estimated by a multivariate version of Jackson theorem 
(see, for example, [15]). which shows that E n (f) = 0(n~ p ) for / £ C P (Q), 
p £ R + . It becomes an effective approximation tool in the uniform norm 
when its operator norm (the so-called Lebesgue constant) grows slowly (cf. 
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[TBI [T51 ITT1 [S])- The hyperinterpolation has been used effectively in several 
cases: originally for the sphere [TH1 EH] , and more recently for the square [H[S], 
the disk [TTJ, and the cube [6]. We will use our new cubature formulae to 
construct a hyperinterpolation operator of three variables for the Chebyshcv 
weight function on the cube. We show that the computation can be carried 
out efficiently using the 3-dimensional FFT and that the algorithm can be com- 
pletely vectorized. We will also present numerical results on hyperinterpolation 
of several test functions. 

The paper is organized as follows. In Section 2 we construct new cubature 
formulae and report results of numerical tests, where comparisons are made 
with tensor-product Gauss-Chebyshev formulae. Hyperinterpolation in three 
variables is considered in Section 3, where we show how to compute it effec- 
tively and report the results of numerical tests. Finally in Section 4, we obtain 
a new (nontensorial) Clenshaw-Curtis type formula in the cube by integrating 
the hyperinterpolant in Section 3 and show that it has a clear superiority over 
tensorial Clenshaw-Curtis and Gauss-Legendre cubature on nonentire test inte- 
grands, a phenomenon known for 1-dimensional and 2-dimensional Clenshaw- 
Curtis formulae (see [2"UIIT9"] ). 

2 Algebraic cubature for the rf-dimensional Cheby- 
shev measure. 

We consider cubature formula for the product Chebyshev weight function Q, 
which is normalized so that its integral over [— 1, l] d is 1. For d = 1, we write 
w(x) — W\(x). 

Let 11^ denote the space of polynomials of total degree < n in d variables. 
We write II„ if d = 1 . The Gaussian quadrature formula for w takes the form 



(9) 



r 1 1 " 

/ f(x)w(x)dx =-Y,f {cos , V/gIWl 

11 k=l 

For d — 2, a cubature formula of degree 2n — 1 needs at least (cf. [13"] ) 

^ = di m( nL l) + [|j=^±ll + [fj (io) 

many nodes. Cubature formulae that attain this lower bound can be constructed 
for the product Chebyshev weight W 2 {x) (see and the references therein) 

by studying common zeros of associated orthogonal polynomials. In [TJ, these 
cubature rules were derived by an elementary method which depends on a fac- 
torization of the Gauss-Lobatto quadrature into two sums, over even indices and 
odd indices, respectively. This factorization method was also used for d > 2 in 
[TJ and yields a cubature formula of degree 2n — 1 for Wd with roug hly n d /2 d / 2 
many nodes. 

A close inspection of the factorization method shows that it actually allows 
us to derive cubature formulae of degree 2n — 1 for Wd with roughly 2(n/2) d 
many nodes. This number of nodes is substantially less than n d of the product 
cubature formula or n d /2 d / 2 of the formulae in [T], although it likely far from 
optimal as seen from 
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We start with the Gauss-Lobatto formula for w on [—1,1]. It takes the form 

/ 1 i /(^)^)^-^2 /(_i)+ £ /(cos ^ )+ ^ /(i) ) (n) 

which again holds for all / € H2 n -i- We proceed to factor this rule into two 
terms. The factorization depends on whether n is even or n is odd. Define 

^/:=M^(-l)+E/( C0S? f)+^( 1 ) 
n = 2m : V i_1 / (12) 



'°/:='£/(cos^) 



n 



and define 



C/'= ■ I >./(-s 2 f) + i/(l) 



1 / m_1 

n = 2m-l: ; ' N (13) 

/ m— 1 

^/:=Mi/(-l) + E/( C0S ^ 

where we use the superscripts £7 and O to signify that the sum is taken over even 
indices or odd indices, respectively. Evidently, the quadrature (11) becomes 

el 

f(x)w(x)dx = i*f + j°f , v/ e n 2 „_! , 

i 

by definition. 

The Chcbyshcv polynomials, T n , are orthogonal with respect to w on [—1,1], 

T n (t) := cosn9 , t = cos 9 . 

The following elementary lemma plays a key role in constructing cubature for- 
mulae on [—1, l] d . 



Lemma 2.1 For n > and k G Z. 



, i 0. k ^ mod n , , i , 

IuT k = { i , n anrf /°T fe = {I, fc = 0,2n,4n 

i , fc = mod n 1 - 



0, fc 7^ mod n 
fc = 0,2n,4n 
|, fc = 0, n, 3n, 



Proof. The proof follows from elementary trigonometric identities. For exam- 
ple, for n = 2m, an elementary calculation shows that 

jo Tk = l _y cosk i2^ sin far - sinfc7r 



n-f-f zm 4msin|^ 2nsin^ 
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from which I^Th = for k ^ mod n follows immediately. The case when k 
is a multiple of n follows from the first equal sign of the above equation without 
summing it up. Similarly, 



1 / 1 



cos kn + cos k^ 



1 \ sin kir cos — 



3=1 



2n sin 



from which the stated result follows. The proof for n = 2m — 1 is similar and is 
omitted for brevity. q.e.d. 

Let a £ {E, 0} d , that is, 

a = (o"i , . . . , <Jd) with <7j = E or <7j = O. 
For a function / : W l i— ► R, we define the sum 

/■"■i ... A CT ti f 

n n J 

as a d-fold multiple sum in which J CTfc is applied to the fc-th variable of /. Let 
us define 

- E a i = 
O a^E 1 ' 

For each a £ {E, 0} d , we then define 

Since the sum introduces a symmetry among a £ {E, 0} d , there are 2 d ~ 1 dis- 
tinct I° d f sums. 

Theorem 2.2 For e? > 1 and each a £ {E, 0} d , the cubature formula 

f(x)W d (x)dx = 2 d - l I°J (15) 

-1,1]* 

is exact for f £ n d n _ 1 and its number of nodes, N, satisfies 



N 



([|J) d (l + o(n- 1 )). 



Proof. For k = (h, . . . , k d ) £ N d let T k {x) := T kl ( Xl ) ■ ■ ■ T kd {x d ), which is a 
polynomial of total degree |fc| :— k\ + ■ ■ ■ + k d . It suffices to establish (15) for 
/ £ {T k : \k\ < 2n — 1}, since this set is an orthogonal basis of LT^. In this case 
we have 



[-1,1? 



T k (x)W d (x)dx = 2 d - x 



flT, . . . J a drp I flT, . . . fdT, 



By the orthogonality of T k , the left hand side is equal to 1 if k = (0, . . . , 0) 
and zero if k ^ (0, . . . , 0). From the definition of and 1°, it is evident that 
1^1 = I°l = 1/2. Hence, for k = (0, ... ,0), the right hand side is equal to 
2 d-i(2-d + 2 -d) = 1; verifying the equation for k = (0, . . . , 0). 
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Assume now < \k\ < 2n — 1. If one of ki ^ mod n, then I^ d Tk = 



by Lemma 2.1 We are left with the case that ki = mod n for all i. Since 
\k\ < 2n — 1, there can be at most one ki = n. Furthermore, |fc| > shows that 
there is exactly one ki = n. Thus the right hand side becomes /^T„ + I%*T n — 
I%T n + I°T n , which is zero as I%T n = 1/2 and I°T n = -1/2 according to the 
Lemma [2~Tj q.e.d. 



For the case of d — 2, Theorem |2.2| contains two distinct cubature formulae 
for a = (E, E), (E, O), respectively, whose number of nodes are either equal to 



N* in (10 1 or N* + 1, those are the ones that have appeared in [HJ|22], and 
later in pQ, as mentioned earlier. For d = 3, there are 4 distinct formulae for 
a = (E,E,E), (E,E,0), (E,0,E), (0,E,E), respectively. For n = 2m, the 



number of nodes is 



for cr = (E, E, E) and 



N 



N = 



(n+ 1) 3 + (n+ 1) 



(n+l) 3 -(n+l) 



for cr= (E,E,0), (E,0,E), (0,E,E), respectively. 

In order to demonstrate the effectiveness of the new cubature formula, we 



present in Figures 1-2 numerical results of (15) with a = [E,E,E) on the inte- 
grals of six text functions with respect to the product Chebyshev measure on the 
3-cube. The first three functions are analytic entire (a polynomial, an exponen- 
tial and a gaussian), whereas the other three are less smooth: one analytic but 
not entire (a 3-dimensional version of the Runge test function) , one C°° nonan- 
alytic, and one C 2 . These functions are analogues of test functions for algebraic 
cubature in dimension 1 and 2, see [201 US]- We compare them with two natural 
choices for cubature on a tensor product domain: the tensor-product Gauss- 
Chebyshev and Gauss-Chebyshev-Lobatto formulae. The results, obtained with 
Matlab (cf. [10]), demonstrate the superiority of the new formula in all cases, 
especially for the less smooth functions, in terms of number of function evalua- 
tions. It should be pointed out that, however, the superiority for the less smooth 
functions arises for even n (a sort of parity phenomenon) . Other numerical tests 
(not reported for brevity) have shown that the cubature formula has the same 
behavior for a =(E,E, O), (E,0,E), (0,E,E). 

A natural question associated with cubature formulae is polynomial interpo- 



lation. Let X n -\ denote the set of the nodes of the cubature formula (151. The 
interpolation problem looks for a polynomial subspace, S , of the lowest degree 
such that 

P(x) = /0), x e X n _ 1; V/ g C(R d ) 

has a unique solution in S. In the case of d = 2, this problem is completely 
solved in [22], where S is a subspace of II 2 which includes and compact 

formulae of the fundamental interpolation polynomials are also given there. For 
d > 2, however, the problem is much harder, since the number of nodes of 
our cubature is far from dim(II^). For example, if d = 3, then dim(II^_ 1 ) = 
n(n + l)(n + 2)/6 ~ n 3 /6, whereas our cubature has « n 3 /4 many nodes. The 
problem essentially comes down to study the polynomial ideal that has X n _i 
as its variety (see [2"3"]1. 



G 



A simpler approach to polynomial approximation via these new nodes is 
given by hyperinterpolation, as described in the Introduction. In the next sec- 
tion we shall apply such a method in the 3-dimensional case. 



3 Implementing hyperinterpolation in the 3-cube. 

We now use cubature formula ( |15| to construct hyperinterpolation as in 

for 

the 3-cube il = [— 1, l] 3 . In this case, {p a } is the product Chebyshev orthonor- 
mal basis (cf. j9]), i.e. 

p a (x) = T ai (x 1 )T a2 (x 2 )T a3 (x 3 ) , (16) 
where Tfc(-) = v2 cos(fc arccos(-)), k > and Tq(-) — 1. Moreover, let 

C n = I cos — , k = 0, n 
I n 

be the set of n + 1 Chebyshev-Lobatto points, and C„, C° its restriction to 
even and odd indices, respectively. Then, 

X n = (Cn+l X X ^n+l) U \Pn+l X ^n+1 X ^n+1 ) I (1?) 



with (cti, (T2, f"3) S {£7, O} 3 , see (14 1. The weights of the cubature formula (15 1 

for £ G -X" n , are 



1 if £ is an interior point 

, 1/2 if £ is a face point 

w £ ~ ( n _)_ ^3 | 1/4 if £ is an edge point 

1/8 if £ is a vertex point 



Note that, since 

dim(Tl 3 (f2)) = (n+ l)(rc + 2)(ra + 3)/6 < N = card(X„) w n 3 /4 , 

the polynomial £ n / in ^ is not interpolant. 
Now, defining 

m) = F(b, = < (19) 

[ £ G (C„+i x C„+i x C„ +1 )\X„ 

we can write 

c « = E w af(0Pa(0 

= E ( E ( E ^(6,6,6)r Ql (6) f Q2 (6) pf Q3 (6) 

6eC„ +1 \6ec„ + i \6£C„ + i / / 



' 3 \ n+1 /n+1 /n+1 i \ ■ \ ■ 

vs=l / »=0 \i=0 \fc=0 



cos - cos 



n + 1 / n+1 
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Cubature rel. errors for exp(x+y+z) 
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Figure 1: Relative cubature errors versus the number of function evaluations 
for three test functions. 
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Cubature rel. errors for (xW+z 2 ) 3 ' 2 



-© — Our formula 

— * — Tens. -prod. Gauss-Chebsyhev 
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1 1.5 2 

Number of function evaluations 



x 10 



Figure 2: Relative cubature errors versus the number of function evaluations 
for three test functions. 
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where 



a = (ai,a 2 ,a 3 ) G {0, 1, . . • ,n} 3 , [3 as = i ^ ^" = ' s = 1 ' 2 ' 3 - 

This shows that the 3-dimensional coefficients array {c a } is a scaled Discrete 
Cosine Tranform of the 3-dimensional array 

/ in jir kit \ 

F ijk = F cos — — -,cos — — -,cos — — , 0<z,j,fc<n, (20) 
\ n + 1 n + 1 n + 1 / 

where we eventually pick up only the (n + l)(n + 2)(n + 3)/6 ~ n 3 /6 hyperin- 
terpolation coefficients corresponding to \a\ = a\ + a2 + a-^ < n. 

A fast implementation of hyperinterpolation is now feasible (for example 
in Matlab), via the FFT. Indeed, we have written a Matlab code (see [5]), 
completely vectorized by several implementation tricks, whose kernel can be 
summarized as follows: 

Algorithm: Fast total degree hyperinterpolation in the 3-cube 

(i) construct the hyperinterpolation point set X n as union of the two subgrids 
in pi; 



(ii) compute the cubature weights in ( |18| ; 

(iii) compute the 3-dimensional array {-Fyfc} at the complete grid C n +i x 



C„+i x C n +i by (19) (notice that / is evaluated only at X„ 



(iv) compute the 3-dimensional array of coefficients {c a } by three nested ap- 
plications of the 1-dimensional Real(FFT(-)) operator; 

(v) select the coefficients {c Q } corresponding to the triples a — (ai, 02,0(3) 
such that \a\ — a\ + a 2 + «3 < n. 

We recall that there is a simple way to approximate a function in the 3-cube 
by tensor-product of polynomials of degree n, that is, by a tensor-product dis- 
crete Chebyshev series (ultimately a tensor-product hyperinterpolant). Such an 
approximation uses (n + l) 3 function evaluations, and (n + l) 3 coefficients. In 
contrast, let us stress again the following facts on our total-degree hyperinter- 
polation of degree n in the 3-cube: 

Remark 

• the number of hyperinterpolation nodes, or function evaluations, is equal 
to card(A„) w n 3 /4; 

• the number of hyperinterpolation coefficients is dim(II 3 ) n 3 /6. 

In order to compare the performances of total-degree and tensor-product 
hyperinterpolation in the 3-cube, we show, in the following figures, the hy- 
perinterpolation errors versus both the number of nodes and the number of 
coefficients on the six test functions already used in Section 2, and we choose 



again (ai , 02 , 03) = (E,E,E), see (17 1. The errors are relative to the maximum 
deviation of the function from its mean and are computed on a uniform control 
grid. Since the computation of the coefficients via the FFT has roughly the 
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same cost for both kinds of hyperinterpolation, we have chosen the number of 
function evaluations as a measure of computational cost for the construction, 
and the number of coefficients as a measure of the compression capability of the 
algorithms. 

The situation here is in some sense opposite to that of Figures 1-2. In- 
deed, total-degree appears superior to tensor-product hyperinterpolation on the 
smoothest functions, but not on the less smooth ones. As it is natural from the 
observation above, the behavior of total-degree hyperinterpolation in terms of 
number of coefficients is better than that in terms of number of nodes (function 
evaluations). 

4 A Clenshaw-Curtis-like formula in the cube. 

In the recent paper [TS], perusing an idea already present in [T7], it has been 
shown how hyperinterpolation allows us to construct new cubature formulae. 



Given h £ L 2 d (fl) and / <E C(f2), we can approximate the integral of hf in d/j, 



as 



h(x) f(x) djiK I h(x) C n f(x) djjL 



= E c « m «= E < ( 21 ) 

\a\<n ZeX n 

where the generalized "orthogonal moments" {m a } and the cubature weights 
{A^} are defined by 

m a := / h(x) p a (x) dfi , Aj := w 5 ^ Pa(£)rn a . (22) 



Observe that the cubature formula (21 1 is exact for every / <E II„(f2), and that 



{m a } are just Fourier coefficients of h with respect to the /i-orthonormal basis 
{p a }- 

Concerning stability and convergence of such cubature formulae, the follow- 
ing result has been proved in |19j : 

Theorem 4.1 Let all the assumptions for the construction of the cubature for- 



mula (21) be satisfied, and in particular let h £ L^„(f2). Then the sum of the 
absolute values of the cubature weights has a finite limit 



E l A ?l= / l*0=)|dM. (23) 



Notice that ( 23 1 ensures that the sum of absolute values of the weights is 



bounded, and thus by recalling that C n is a projection operator on 11^ (f2) we 
obtain the Polya-Steklov type (cf. [H]) convergence estimate 



h{x) f{x) d[i 



E 



< 



\h(x) \ d/i 



sup V 

n ^x n 



(24) 
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Hyperinterpolation rel. errors for (x+y+z) 
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Figure 3: Hyperinterpolation relative errors versus the number of function eval- 
uations for three entire test functions. 
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Hyperinterpolation re!, errors for 1/(1+16(x 2 +y 2 +z 2 )) 
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Figure 4: Hyperinterpolation relative errors versus the number of function eval- 
uations for three noncntirc test functions. 
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Hyperinterpolation rel. errors for (x+y+z) 




-Tens. Prod. Hyperinterp 
Tot. Deg. Hyperinterp. 



yi o _ 
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Number of hyperinterpolation coefficients C a x 10 4 



Hyperinterpolation rel. errors for exp(x+y+z) 
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Figure 5: Hyperinterpolation relative errors versus the number of hyperinterpo- 
lation coefficients for three entire test functions. 
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Hyperinterpolation re!, errors for 1/(1+16(x 2 +y 2 +z 2 )) 
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Hyperinterpolation rel. errors for exp(-1/(x 2 +y 2 +z 2 )) 
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Figure 6: Hyperinterpolation relative errors versus the number of hyperinterpo- 
lation coefficients for three noncntire test functions. 
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where E n (f) denotes the error of the best polynomial approximation of degree 
n to / in the uniform norm. 



Now, applying ( 21 1-( 22 ) in the case 



d[i = w(x)dx, w € LjJCl) , with h = — € L l d {Vt) , (25) 

w 

(since then h 2 = l/w 2 e L^Q)) we obtain, via hyperinterpolation, a cubature 
formula for the standard Lebesgue measure from an algebraic cubature formula 
for another measure (absolutely continuos with respect to the former). The 
specialization of this approach to the 1-dimensional Chebyshev measure gives 
ultimately the popular Clenshaw-Curtis quadrature formula [7 . An extension 
to dimension 2 has been studied in |19j . Here we apply the method in dimension 
3, obtaining a new nontensorial Clenshaw-Curtis-like cubature formula in the 
3-cube. 

In Figures 7-8 we display the relative errors of such a formula for (ci , 02 , 03 ) = 



(E,E,E) (cf. ( IT I ) on the six test functions already used above, compared 
with those of the tensor-product Clenshaw-Curtis, Gauss-Legendre, and Gauss- 
Legendre-Lobatto formulae. The numerical results have been obtained with 
Matlab, using jTUj for the Gaussian formulae and [2T] for the tensor-product 
Clenshaw-Curtis formula. 

In particular, we see that with the entire test functions nontensorial Clenshaw- 
Curtis cubature is more accurate than the tensor-product version, but less ac- 
curate than the other two tensor-product formulae. On the other hand, in the 
less smooth cases the nontensorial Clenshaw-Curtis formula is better than all 
the other three, especially for odd hyperinterpolation degrees n, which corre- 



spond to use n+l even in ( 15 1 (again a sort of parity phenomenon, cf. Figure 
2). This behavior echos that of 1-dimensional and 2-dimensional Clenshaw- 
Curtis formulae (see [2(3 [19]). Other numerical tests (not reported for brevity) 
have shown that the other versions of the nontensorial Clenshaw-Curtis formula, 



corresponding to (ci, cr 2 , 03) = (E,E,0), (E,0,E), (0,E,E) in (17), produce 
essentially the same results. 
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